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Abstract 

In recent years, many types of elliptical Radon transforms that integrate functions over 
various sets of ellipses/ellipsoids have been considered, relating to studies in bistatic synthetic 
aperture radar, ultrasound reflection tomography, radio tomography, and migration imaging. In 
this article, we consider the transform that integrates a given function in R ra over a set of ellipses 
(when n = 2) or ellipsoids of rotation (when n > 3) with foci restricted to a hyperplane. We 
show a relation between this elliptical Radon transform and the regular Radon transform, and 
provide the inversion formula for the elliptical Radon transform using this relation. Numerical 
simulations are performed to demonstrate the suggested algorithms in two dimensions, and these 
simulations are also provided in this article. 
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1 Introduction 


Radon-type transforms that assign to a given function its integrals over various sets of ellipses/ellipsoids 
arise in migration imaging under an assumption that the medium is acoustic and homogeneous. The 
aim of migration is to construct an image of the inside of the earth from seismic reflections recorded 
at its surface [21 [25]. A graphical approach called classical migration was developed systematically 
by Hagedoorn mi. Classical migration had been abandoned after the wave-equation method was 
introduced by Claerbout J5] in 1971. Gazdag and Sguazzero pointed out that the classical migration 
procedures that existed at that time were not based on a completely sound theory [9] . However, 
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a correct construction for the wave-equation method was often difficult to find because the ex¬ 
periment data did not fit into a single wave equation. To adapt the diffraction stack to borehole 
seismic experiments, a new approach to seismic migration was found. This approach gave classical 
migration a sound theory. After that discovery, classical migration has attracted many researchers 
in the field. The underlying idea is that seismic data in the far field can be regarded as if the 
data are coming from integrals of the earth’s acoustic scattering potential over surfaces determined 
by the velocity model M- These Radon-type transforms relate to migration imaging as well as 
Bistatic Synthetic Aperture Radar (BiSAR) [2j [6j s f3j[T4i[32], Ultrasound Reflection Tomography 
(URT) |ll m 26], and radio tomography [271 2h- 29] - 

Because of these applications, there have been several papers devoted to the topic of elliptical 
Radon transforms. The family of ellipses with one focus fixed at the origin and the other one 
moving along a given line was considered in HU. In the same paper, the family of ellipses with 
a fixed focal distance was also studied. The authors of mm dealt with the case of circular 
acquisition, when the two foci of ellipses with a given focal distance are located on a given circle. 
A family of ellipses with two moving foci was also handled in [6]. Radio tomography, which uses a 
wireless network of radio transmitters and receivers to image the distribution of attenuation within 
the network, was discussed in [2711281 [29]. They approximated the obtained signal by the volume 
integral of the attenuation over this ellipsoid. One work H7UH] derived two inversion formulas of 
this volume integral of the attenuation over this ellipsoid and studied its properties. Many works 
found an approximate inversion for elliptical Radon transforms. 

Here we consider migration imaging and introduce a new type of an elliptical Radon transform 
obtained by restricting the position of the source and receiver in migration imaging. We find an 
explicit inversion formula for this elliptical Radon transform arising in migration imaging which is 
the line/area integral of the function over the ellipse/ellipsoid with foci restricted to a hyperplane. 

The rest of this paper is organized as follows. The problem of interest is stated precisely and the 
elliptical Radon transform is formulated in section [2] In section 3, we show how to reduce the ellip¬ 
tical Radon transform to the regular Radon transform. The numerical simulation to demonstrate 
the suggested 2-dimensional algorithm is presented in section 4. 


2 Formulation of the problem 


Let s£R 3 and r € R 3 represent 3-dimensional source and receiver positions, respectively. For fixed 
points (s,r), an isochron surface I( s ,r,t) 0 is a surface consisting of image points x = (xi,X2,X3) £ R 3 
associated by the travel time function r(x, y), which gives the travel time t £ [0,oo) between two 
points x £ R 3 and y £ R 3 with a known velocity u(x) (see BED- 

Mathematically, I( s ,r,t) can be described as the set of the image points x satisfying the constraint 
that the total travel time from the source s though the image point x to the receiver r is constant 

1 Hagedoorn called this set a surface of maximum concavity (see ini '). 
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and equal to t. The isochron surface 

^(s, r ,i) = 


I(a,r,t) can be represented as 
{x £ M 3 : t = t(x, s) + r(x, r)}. 


Seismic experiments yield data g{ s, r, t) which are functions of the source position s, receiver posi¬ 
tion r, and time t. Assuming an object function /(x) on R 3 , the data g is modelled by the integral 
of / over i-e., 

g(s,r,t) = j f (x)dx (see UB]). 

^(s,r,t) 

The two detectors have nonzero sizes and time is also passing, so it is reasonable to assume that 
what we measure is the “average”-concentrated near the location of the detectors and nearly zero 
sufficiently far away from the detectors-over a small region of space and a small time interval 
preceding the time t. In mathematical terms, our data can be written as 

g(s,r,t)= / /(x)d(r(x, s) + r(x, r) — f)dx, 

]R 3 

where 5 is the Dirac delta function. (Most works [II l 2 l l 6 l [TOl [T3l HH] dealing with an elliptical Radon 
transform consider the arc length measure instead.) Miller, Oristaglio, and Beylkin suggested this 
model and approximately recovered the object function / by appropriately weighted back projection 
of the data in [16]. However, we present the explicit formula for reconstructing / using the specific 
minimum data of g by changing the variables of positions of the two detectors. 


The total dimension of the data g is 7. To reduce the overdeterminacy, we assume that s and r 
are located on a hyperplane, say X 3 = 0 and also that the difference vector s — r is parallel to a 
given line, say the aq-axis. (In 2 dimensions, the difference vector s — r is automatically parallel to 
the aq-axis because the hyperplane is the line.) Hence, for any u 2 £ R, we write s = ( 5 ,^ 2 ,0) and 
r = (r, ?/ 2 , 0 ). Also, assuming a velocity u(x) is set to a constant 1 , r(x, y) becomes the distance 
|x — y| between x and y. Since x £ I(s,r,t) satisfies t = |x — (s, 112 , 0)| + |x — (r, tt 2 , 0 )|, a point 
x € I( t) can be described by the ellipsoid equation 

(aq - (s + r)/2) 2 (x 2 - u 2 ) 2 x% = 

f 2 /4 t 2 /4 - (r - s) 2 /4 t 2 /4 - (r - s ) 2 /4 

Let us choose a constant a between 0 and 1. We consider the fixed value of r satisfying the equation 
r = s + to, which is a new additional requirement. Then a point x £ 7( s , r ,i) can be also described 
by the equation 

(aq - (s + r)/ 2) 2 (x 2 - u 2 ) 2 xj = 2 

1/4 (1 - a 2 )/4 (1 - a 2 )/4 

If we set u\ = (s + r)/ 2, then our data become 


g((ui - ta, it 2 , 0), (ui + ta, u 2 ,0),t) 


= / /<x)s (/ 

R3 X 

:= R Et af(u,t), 


(xi - til ) 2 (x 2 - U 2 ) 2 

1 / O \ / J ' 


1/4 


(1 — a 2 )/4 (1 — a 2 )/4 


— t dx 
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where u = (ui,U 2 ) € M 2 . Here, we restrict the positions of the source and receiver in migration 
imaging, say, (u\ — ta,U2,0) and (u\ + ta,U2,0) for u £ JR 2 . In general, reducing the space of 
positions for two devices is more useful and economical. If the function is odd with respect to 
X 3 , then R.E,af is equal to zero. We thus assume the function / to be even with respect to X 3 : 
/(x 7 ,^) = /(x 7 ,—x 3 ) where x = (x 7 ,X 3 ) £ l 2 x 1. We call RE,af an elliptical Radon transform, 
since this is the surface area integral of / over the set of these ellipsoids. We generalize this 
3-dimension setup to n-dimension and define a more general form of an elliptical Radon transform. 

Definition 2.1. Let ai, 02 , • • • ,a n >0be given numbers, let A := diag(a \, • • • , a n ) denote the n x n 
diagonal matrix with diagonal entries ai, and let f be a locally integrable function on M n , even with 
respect to x n . The elliptical Radon transforms of f is defined by 

RE,Af(u,t ) = J f(x)S (|R _i ( x - (u, 0 ))| - t) dx, for (u ,t) £ x [ 0 ,oo). 

R n 


Note that we do not need the condition a* < 1 any more. If all ai, i = 1,2,--- ,n, are equal to 1, 
the elliptical Radon transform Re,a} becomes the spherical Radon transform with the centers of 
the sphere of integration located on the hyperplane, a well-studied transform (see [312 0 El m 

□21DEU2Q1 EH Eg EHUMl EH EE] ). 

3 Inversion formula 


Here we assume that the object function / does not touch the detectors; that is, the support of / 
does not intersect the hyperplane x n = 0 where the two detectors are located. 


To obtain an inversion formula for the elliptical Radon transform, we manipulate Re,A$■ By 
changing the variables H“ 1 (x — (u, 0 )) —> x, we can write 


R E ,Af ( u , t) = |a|! J f(Ax + (u, 0))d(|x| - f)dx 

R n 

= |a|R” -1 f f (Ayt + (u, 0))dS(y) 

|y|=i t 

= 2|a|it" -1 J f(A(y', y 7 1 - [y 7 | 2 )f + (u, 0)) 


(1) 


|y'l<i 

where y = (y 7 , y n ) £ M n and |a|i = aqa 2 • • • a n . Here dS( y) is the area measure of the unit sphere. 


Let 2t = {z = (z ',z n ) = (zi,Z 2 ,--- ,z n -i,z n ) € M n : 0 < |z 7 | 2 < z n j and 53 = {x = (x 7 ,x n ) = 
(xi,X 2 , ■ ■ ■ , x n -i,x n ) € M n : x n > 0}. We define a map m n) A : 2t 5S by 

m n}A { z) = (• Az',a n sjz n - |z 7 | 2 ), 

where A := diag(a\, ■ ■ ■ a n - 1 ) is the n-lxn -1 diagonal matrix. 
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Proposition 3.1. 


The map m n ^A '■ 21 —> 53 is a bijection with the inverse map m nj A 1 : 55 —>• 21 defined by 

m B ,A _1 ( x ) = (A _i x 7 , |-4 _i x| 2 ). 


We /iawe 


{m n ,A 1 (x)eM n :[ J 4 x (x - (u, 0))| = t, x £ 53} 
(-2i" 1 u,l) t 2 — j^4- 1 u| 2 


= < z £ 21 : z 


Rl( u ) 


Rl( u ) 


w/iere z^(u) = y/l + 4|^4 _1 u| 2 . 


The map rri n _A 1 transforms an ellipsoid into a hyperplane. Changing variables using this map m n _A 
plays a critical role in reducing the elliptical Radon transform to the regular Radon transform. 

Proof. We can easily check that m ni A -1 o m njJ 4 (z) = z for z £ 21 and m n ^ o m n> A _1 (x) = x for 
x £ 23, so m Uj A ■ A —>• IB is a bijection. 

Consider m niJ 4 _1 (x) • (—2R _1 u, lj/z/^u) for x £ 53 with A _1 (x — (u, 0))| = t: 

m n ,A _1 (x) ■ = (A-V, |^4- 1 x| 2 ) ■ 1U ’ ^ 


_ Ri( u ) 

—2v4” 1 x' • R _r u + |R _1 x| 2 

^(u) 

|R _1 (x - (u, 0))| 2 - |A _1 u| 2 _ t 2 - IR-iu' 2 


v d( u ) 


Rl( u ) 


ri( u ) 


( 2 ) 


□ 


We define the function k{ z) on M n by 


f °m n , A { z) 


if 0 < |z | 2 < z n , 


fc ( z ) = ^ V-n - |z'| 2 

0 otherwise, 

where z = (z 7 , z n ) £ M n . This is equivalent for > 0, to 

/(x) = x n a~ l k o m ni A _1 (x) = IW^x] 2 ), 

where x = (x 7 , x n ) £ M n . By the evenness of /, we have 

/(x) = |^n|a^ 1 fe(A _1 x 7 , |A _1 x| 2 ). 


(3) 
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Let the regular Radon transform Rk(eg,s ) be defined by 


Rk(eg, s) = J k(see + rj)drj, 




(eg, s') € S 71 - 1 xR, 


where s £ R and for 6 = ( 8 \, 82 , ■ ■ ■ , 0 n -i) £ [0,27r) x [ 0 , ir\ 


n —2 


eg = 


/ sin 9 \ sin #2 • • • sin 0 n _i \ 
cos 9 1 sin $2 ■ ■ ■ sin 9 n -\ 
cos 62 sin #3 • • • sin 9 n -i 

cos Q n - 2 sin 0 n _ 1 
cos 0 n _l 


e 


n—1 


This can be represented by 

Rk(eg, s ) = 


1 


COS 0 n _l| 


k [ z', - 


en • z' 


+ 


COS0 n _i COS0 n _i 


dz 7 , 


(4) 


for 0 n _i 7 ^ 7t/2 and eg = (eg,cos0 n _i) € <S n 1 . The following theorem shows the relation between 
R e ,aJ and Rk. 


Theorem 3.2. Let f e C(M n ) 6 e even in x n and have compact support in M n 1 
we have 


R E ,Af( u,t) 


2 | a | i ^ / (—2x4 1 u, 1) t 2 - 1 u| 2 \ 

rj( u ) V rj( u ) ’ *Li( u ) / 


xl\{ 0 }. 


Then 


(5) 


Again, 1 / 4 ( 11 ) = \/l + 4| J 4~ 1 u| 2 . 


Proof. Combining two equations (pQ) and ([3]) , we have 

R E ,Af(n,t) = 2 t n |a|i J k(y't + A _ 1 u, |y'f + A _ 1 u | 2 + (1 - |y'| 2 )t 2 )dy'. 
|y'l<i 

Let us consider the variable change 


z' = y't + T _1 u 




The Jacobian of this transformation is t 1 n , so 


R E}A f(n,t ) = 2|a|if 


k( 1 


Iz'l 2 + 1 2 


| z'— A~ 


zl — A 1 u| 2 )dz / . 
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Since k has compact support in {z £ l 2 : |z 7 | 2 < z n }, we have 

Re,aI (u, t) = 2|a|it J k(z', |z 7 | 2 + t 2 — |z 7 — T _ 1 u| 2 )a!z 7 

R ™- 1 

= 2|a|it J k{z', t 2 — |j4 _1 u| 2 + 2 (yi~ 1 u) • z')dz'. 
i "- 1 

We recognize the right hand side as the integral along the hyperplane perpendicular to 

(- 2 A“ 1 u,l)/z/ yi (u) 

with (signed) distance 

0 1 2 - \A~ 1 u\ 2 )/v a (u) 

from the origin. In this case, the measure for the hyperplane becomes z/|(u )dz'. Setting 
eg = (— 2A~ 1 u 1 l)/z^(u), and s = (t 2 — |^4 _ 1 u| 2 )/^(u) 
in equation Q we have the desired formula from equation ([ 6 ]) . 


( 6 ) 


□ 


Using the projection slice theorem for the regular Radon transform, we obtain an analog of the 
projection slice theorem for the elliptical Radon transform Re,a^'- 

Theorem 3.3. Let f £ C'°°(R n ) be even in x n and have compact support in M n_1 xl \ {0}. Then 
we have for (a, (3) £ M n_1 x M, 

OO _ 

k(oc,/3) = j R E A f e - ^* 2 dt 

o 

where k is the n-dimensional Fourier transform of k. 


Proof. The projection slice theorem implies 

k(aeg) = J Rk(eg, s)e~ lS(T ds. 

R 

_^t 2 — \A~ *u | 2 ^ 2 | —_i 2 

To use this theorem, we multiply equation m by e a ( u) and integrate —from 
— \A~ l u\ 2 /v a (\T) to infinity: 

J Rk( { ~ 2A luA \ t2 ~\ A [ 

J V U4l u ) ^a( u ) 


U 


■ t 2 - |A~ 1 u|“ 
1 " A ( u ) 


t 2 - IA-^i 
"a( u ) 
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Since k has compact support in {z G W 1 : \ 
distance from the origin less than 
/c(cr/^( u )(—2 t4 _ 1 u, 1 )) is equal to 


z'| 2 < z n }, the plane perpendicular to ^ ^(u) 1 ’ 1 ^ w hh 
does not intersect the compact support of k. Hence, 


Rk 


-lA-iuP/i/^u) 




v . d( u ) ’ "a( u ) 


/ 


u A (u) „ „ , - n 2 -i^ u ' 2 CT , A 2 - IH-iul 2 

A -R E ,Af( M)e ^ (u) cC 1 1 


2 |a|it 


-lA-iuP/^Cu) 


F4( u ) 


1 

l a li 


R E ,Af(u,t)e 1 1 W u) 


- 

^a( u ) / 


where in the second line, we used Theorem 13.21 Setting a = —2cjH 1 11 / 1 / 4 ( 11 ) G M n 1 and /3 = 
< 7 / 1/4 (u) G M completes our proof. □ 


Taking the inverse Fourier transform of a function k and using equation Q, we have the following 
inversion formula: 


Corollary 3.4. Let f G C'°°(R n ) be even in x n and have compact support in M n 1 x R\{0}. Then 
/(x) can be reconstructed as follows: 


(-ir-Vnj f 
27T n |a| 2 J 

R n_1 


A t{R E ,Af){0L, t)\ t =\A-i-{(x-x.')\ 2 +a- 2 xl d(X ' 


(7) 


where 

OO 

A Mt) = J J e^-^lP^h^drdp. 

R 0 


Proof. We get for x n > 0, 


/( X ) I 2 ) 


Xr, 


(2vr ) n a r 



0 iA 1 x'•& i/3\A 1 x|-‘ 


4 Ml 


k(a, (3)d/3dot 




4a 


(27r) n a n |a|i 

R™- 1 R 

where in the first line, we used equation ([3]) and in the last line, we used Theorem 13.31 and 


G(u,w) = R E , A f(u,t)e lwt2 dt. 


o 















The evenness of / in x n gives us 


/(x) 


(27r) n a n |a|i 


4/3 


e ia '( A 1 x.') e ip\A x x | 2 q 


d/3da ' 


Changing the variables —Aa/(2/3) -A- a, /(x) is equal to 


(— 2 ) n ” 1 |x- 



(2-7r) n |a| 2 

R " -1 K 

By definition of G , /(x) can be determined through 

OO 

(-l) n - l \x n \ 


g i/9(|A a| —2(A <x)-(A x?)+\A x| )q (a, fi) iP^dfidi 


a. 


27r n |a| 2 


JPUA-^-A-G'f+a-Gl-t 2 ) 


RE,Af( a it)\/3\ n 1 dtd/3da. 


l n ~ 1 K 0 


□ 


When n = 2, equation 0 becomes 

M 


/(x) = 


27T 2 ( 


i^JJ y" 1 («-^ l) l 2 +« 2 2 a; i- t 2 ) J R jE ;^/(Q ; , 


1^21 

7T 2 ofa2 

since in the distribution sense, 


Jfw 


RE,Af( a R) 


(a — x\)\ 2 + a 2 z x 2 — t 2 ) 2 


-dtda, 


j e ixfi 


4 2-dimensional numerical implementation 


Here we discuss the results of 2-dimensional numerical implementations. 
In 2-dimension, equation Q becomes 


Re,aI{ u i t) = 


2aid2t 


— 1„,|2 


,Rk 


{—2a 1 l u, 1) t 2 — \a 1 L u 


y^l + 4|aj~ 1 u| 2 \ ^/l + 4|a]~ 1 tt| 2 * ^/l + 4| 


a, 1 u |2 


Setting (cos 0, sin 0) = (—2a x 1 u, \)/\J\ + 4|a x x u| 2 € S' 1 and 


s = 


t 2 — | a x 1 u | 2 
•^/l + 4|a^ 1 u| 2 


( 8 ) 
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in equation l]8|). we have 


a lC ot0 

—2 ,Ks,0) 


if s esc 6 > 


cot 2 6 
4 


2Rk(cos 0 , sin 0, s) 


esc 6 1 


aid2h(s , 0 ) 


^e,a/ ( - 




where 

/i(s,0) = 

(While /?,/c(sin6 , cos 0, s) is used in Theorem 13.21 Rk(cos 6, sin#, s) is used to utilize the built-in 
function in Matlab. Hence a small change is required.) Again, since k has compact support in 
{z £ M 2 : \z\\ 2 < Z 2 }, the plane perpendicular to eg = (cos 0, sin 9) with distance from the origin 
less than — 4 °s C g does not intersect the compact support of k. Therefore, we have 


s esc 0 + 


cot 2 0 


2 Rk(eg, s ) 


I esc 6» I / aicot6» 

SISSM) - 


if s esc 0 > 
otherwise, 


cot 2 0 
4 


(9) 


First of all, a± and 02 were set to be 0.8 and 1, respectively. In the experiments presented here we use 
the phantom shown in Figure Q] (a). The phantom, supported within the rectangle [—1,1] x [—1,1], 
is the sum of eight characteristic functions of disks. Notice that it has to be even with respect 
to the aq-axis and there are four characteristic functions of disks centered at (0.2, 0.4), (0,0.5), 
(—0.3,0.3), and (—0.5,0.2) with radii 0.2 0.15 0.05 and 0.05, whose values are 1, 0.5, 1.5, and 2 
above the oq-axis. Hence we also include their reflection below the aq-axis. (Actually, our phantom 
has support in 

{x £ M 2 : (aq/oq) 2 + ((aq/aq) 2 + x\) 2 < 1}. 

This implies that k has support in the unit ball and this makes it sufficient to consider the range 
[— 1 , 1 ] in s .) The 256 x 256 images are used in Figured] To reconstruct the image in Figured] 
(b), we have 256 x 256 projections for 0 £ [0, 27t] and s £ [—1,1] in equation ([9]). All projections 
are computed by numerical integration. After finding the function k using the usual inversion code 
for the regular Radon transform, we obtain the function / using equation ([3]). (When using the 
inversion code for the regular Radon transform, the built-in function “iradon” in MATLAB was 
used. The function “iradon” is the inversion of the built-in function “radon” in Matlab which 
considers the number of the pixels where the line passes through. Thus when computing Re,aI , 
we also considers the number of the pixels where the ellipse passes through. We used the default 
version of the function in which the filter, whose aim is to deemphasize high frequencies, is set to 
the Ram-Lak filter and the interpolation is set to be linear.) 


While Figured] (b) demonstrates the image reconstructed from the exact data, Figured] (c) shows the 
absolute value of the reconstruction from noisy data. The noise is modeled by normally distributed 
random numbers and this is scaled so that its norm was equal to 5% of the norm of the exact 
data. In Figure dl (c) the noisy data is modeled by adding the noise values scaled to 5% of the 
norm of the exact data to the exact data. In Figure [2] the surface plots of Figure [1] (a) and (b) 
are provided. Another phantom and reconstruction are shown in Figure [3] The reconstructed two 
dimensional data sets consisting of 256x256 projections using the implemented inversion formula 
have been computed in less than one second (around 0.3 second) on a Intel(R) Cor(TM) i5-3470 
CPU @ 3.20 GHz. 
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(a) 


(b) 


(c) 


Figure 1: Reconstructions in 2 dimensions: (a) the phantom, (b) the reconstruction form exact 
data, and (c) the reconstruction from noisy data 



2-1 



(a) 


(b) 


Figure 2: Surface plots (a) the phantom and (b) the reconstruction from exact data 

5 Conclusion 


This paper is devoted to the study of the elliptical Radon transform arising in seismic imaging. We 
provide an inversion formula for the elliptical Radon transform by reducing this transform to the 
regular Radon transform. Also, we demonstrate our algorithm by providing numerical simulations. 
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Figure 3: Reconstruction in 2D: (a) the phantom and (b) the reconstruction from exact data 
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